#include "StdAfx.h"
#include "RidgeDetector.h"

RidgeDetector::RidgeDetector(void){
	cout << "RidgeDetector initiated without data " << endl;
}

RidgeDetector::RidgeDetector(ImportData* d){
	cout << "RidgeDetector initiated with data: " <<  d << endl;
	this->data = d;
	this->sMin = 0.01; //criteria3: minimum height for a valid ridge
}

void RidgeDetector::setMostConvex(bool c){
	this->m.setConvex(c);
	this->grid_lines2d = this->m.constructRidges(this->grid_points2d);
}

void RidgeDetector::extractRidges3d(){
	cout << "Initiating Ridge Extraction of a 3D Volume.. " << endl;
	Convoluter c(this->data);
	c.initialize3d();
	this->m = Marching(this->data);
	std::vector<MatrixXd> outDx,outDy,outDz,outDxDx,outDyDy,outDzDz,outDxDy,outDxDz,outDyDz;

	//Filter data
	c.computeDerivatives3d(outDx,outDy,outDz,outDxDx,outDyDy,outDzDz,outDxDy,outDxDz,outDyDz);
	//Add found values To DataPoints
	Vector3d g;
	Matrix3d h;
	int x,y,z;
	for(x = 0; x < this->data->getWidth(); ++x){
        for(y = 0; y < this->data->getHeight(); ++y){
			for(z = 0; z < this->data->getDepth();++z){
				if(!this->data->getVolume()[z][y][x].isEmpty()){
					//Map The Found Gradient To Its Valid Point In The Data
					g << outDx[z](x,y),outDy[z](x,y),outDz[z](x,y);
					this->data->getVolume()[z][y][x].setGradient(g);

					//Map The Found Hessian To Its Valid Point In The Data
					h << outDxDx[z](x,y),outDxDy[z](x,y),outDxDz[z](x,y),
						 outDxDy[z](x,y),outDyDy[z](x,y),outDyDz[z](x,y),
						 outDxDz[z](x,y),outDyDz[z](x,y),outDzDz[z](x,y);
					this->data->getVolume()[z][y][x].setHessian(h);
				}
			}//z
		}//y
	}//x

	//Eigen Computation - finding transverse direction and smallest eigenvalue (sigma)
	Matrix3d currentH;
	Vector3d currentG;
	int count = 0;
	double v1_result;
	for(x = 0; x < this->data->getWidth(); ++x){
        for(y = 0; y < this->data->getHeight(); ++y){
			for(z = 0; z < this->data->getDepth();++z){
				if(!this->data->getVolume()[z][y][x].isEmpty()){
					currentH = this->data->getVolume()[z][y][x].getHessian();
					currentG = this->data->getVolume()[z][y][x].getGradient();
					/*cout << "currentG:\n" << currentG << endl;
					cout << "currentH:\n" << currentH << endl;*/

					SelfAdjointEigenSolver<Matrix3d> eigensolver(currentH);
				/*	cout << "The eigenvalues of H are: " << eigensolver.eigenvalues().transpose() << endl;
					cout << endl;
					cout << "The eigenvectors of H are: " << eigensolver.eigenvectors() << endl;
					cout << endl;*/
					if (eigensolver.info() != Success) abort();
						double v1_min; //smallest eigenvalue
						Vector3d v1; //eigenvector corresponding to smallest eigenvalue
						if(eigensolver.eigenvalues()[0] < eigensolver.eigenvalues()[1]){
							if(eigensolver.eigenvalues()[0] < eigensolver.eigenvalues()[2]){
								v1 = eigensolver.eigenvectors().col(0);
								v1_min = eigensolver.eigenvalues()[0];
							}
							else{
								v1 = eigensolver.eigenvectors().col(2);
								v1_min = eigensolver.eigenvalues()[2];
							}
						}
						else{
							if(eigensolver.eigenvalues()[1] < eigensolver.eigenvalues()[2]){
								v1 = eigensolver.eigenvectors().col(1);
								v1_min = eigensolver.eigenvalues()[1];
							}
							else{
								v1 = eigensolver.eigenvectors().col(2);
								v1_min = eigensolver.eigenvalues()[2];
							}
						}
						v1_result = v1.dot(currentG);
						this->data->getVolume()[z][y][x].setEigen(v1); //Assign The  smallest eigenvector to the current DataPoint
						this->data->getVolume()[z][y][x].setSigma(v1_result); //Assign The dot product of smallest eigenvector and gradient as Sigma to the current DataPoint
						//cout << "Result: " << result << endl;
						//double range = 0; //range around 0 that should also be considered a valid ridge
						//if(result >= (0-range) && result <= (0+range)){ //criteria1 eigenvector corresponding to smallest eigenvalue 
						//if(v1_result == 0){ //criteria1 eigenvector corresponding to smallest eigenvalues 
						//	if(v1_min < 0 && v2_min < 0){ //criteria2
						//		//if(this->data->getVolume()[z][y][x].getFTLE() > sMin){ //criteria 3
						//		//	if(!this->data->getVolume()[z][y][x].isRidge()){ //If it is not currently marked as a ridgepoint, mark it.
						//		//		this->data->getVolume()[z][y][x].setRidge(true);
						//		//		count++;
						//		//	}
						//		//}
						//	}
						//}//if result
				}//if !isempty
			} //forz
		} //fory
	}//forx
	cout << "found " << count << " ridgepoints." << endl;
	//March through volume to find ridgepoints
	this->grid_points3d = this->m.march3d();
	//this->grid_lines3d = this->m.constructRidges3d(this->grid_points);
	this->grid_lines3d = this->m.constructSurface(this->grid_points3d);
	//this->m.constructRidges3d(this->grid_points3d);
	//int z;
	//cout << " outDxDy: \n" << endl;
	/*for(z = 0; z < this->data->getDepth(); ++z){
		cout << "Z: " << z << "\n" << endl;
		cout << outDxDy[z] << endl;
		cout << endl;
	}*/
	/*cout << " outDxDz: \n" << endl;
	for(z = 0; z < this->data->getDepth(); ++z){
		cout << "Z: " << z << "\n" << endl;
		cout << outDxDz[z] << endl;
		cout << endl;
	}
	cout << " outDyDz: \n" << endl;
	for(z = 0; z < this->data->getDepth(); ++z){
		cout << "Z: " << z << "\n" << endl;
		cout << outDyDz[z] << endl;
		cout << endl;
	}*/


}

int RidgeDetector::getVertCount(){
	return this->m.vertCount;
}
int RidgeDetector::getTriCount(){
	return this->m.triCount;
}

void RidgeDetector::extractRidges(int slice){
	cout << "Initiating Ridge Extraction of a 2D Slice.. " << endl;
	this->slice = slice;
	Convoluter c(this->data);
	c.initialize2d();
	//Marching march(this->data);
	this->m = Marching(this->data);
	//Marching* march
	//Output Data Collections
	MatrixXd outDx;
	MatrixXd outDy;
	MatrixXd outDxDx;
	MatrixXd outDxDy;
	MatrixXd outDyDy;
	//Convolute The Data
	c.computeDerivatives(this->slice,outDx,outDy,outDxDx,outDxDy,outDyDy);

	//Add Data To DataPoints
	Vector2d g;
	Matrix2d h;
	int x,y;
	for(x = 0; x < this->data->getWidth(); ++x){
        for(y = 0; y < this->data->getHeight(); ++y){
			if(!this->data->getVolume()[this->slice][y][x].isEmpty()){
				//Map The Found Gradients To Its Valid Point In The Data
				g << outDx(x,y),outDy(x,y);
				this->data->getVolume()[this->slice][y][x].setGradient(g);

				//Map The Found Hessian To Its Valid Point In The Data
				h << outDxDx(x,y),outDxDy(x,y),outDxDy(x,y),outDyDy(x,y);
				this->data->getVolume()[this->slice][y][x].setHessian(h);
			}
		}
	}

	//Eigen analysis
	Matrix2d currentH;
	Vector2d currentG;
	int count = 0;
	double result;
	for(x = 0; x < this->data->getWidth(); ++x){
        for(y = 0; y < this->data->getHeight(); ++y){
			if(!this->data->getVolume()[this->slice][y][x].isEmpty()){
				currentH = this->data->getVolume()[this->slice][y][x].getHessian();
				currentG = this->data->getVolume()[this->slice][y][x].getGradient();
				/*cout << "currentG:\n" << currentG << endl;
				cout << "currentH:\n" << currentH << endl;*/

				SelfAdjointEigenSolver<Matrix2d> eigensolver(currentH);
				if (eigensolver.info() != Success) abort();
					double e1 = eigensolver.eigenvalues()[0];
					double e2 = eigensolver.eigenvalues()[1];
					double minEv; //smalles eigenvalue
					Vector2d v; //eigenvector corresponding to smallest eigenvalue
					if(e1 < e2){
						v = eigensolver.eigenvectors().col(0);
						minEv = e1;
					}
					else{
						v = eigensolver.eigenvectors().col(1);
						minEv = e2;
					}
					result = v.dot(currentG);
					this->data->getVolume()[this->slice][y][x].setEigen(v); //Assign The  smallest eigenvector to the current DataPoint
					this->data->getVolume()[this->slice][y][x].setSigma(result); //Assign The dot product of smallest eigenvector and gradient as Sigma to the current DataPoint
					//cout << "Result: " << result << endl;
					//double range = 0; //range around 0 that should also be considered a valid ridge
					////if(result >= (0-range) && result <= (0+range)){ //criteria1 eigenvector corresponding to smallest eigenvalue 
					//if(result >= 0 && result <= (0+range)){ //criteria1 eigenvector corresponding to smallest eigenvalue 
					//	if(minEv < 0){ //criteria2
					//		if(this->data->getVolume()[this->slice][y][x].getFTLE() > sMin){ //criteria 3
					//			if(!this->data->getVolume()[this->slice][y][x].isRidge()){ //If it is not currently marked as a ridgepoint, mark it.
					//				this->data->getVolume()[this->slice][y][x].setRidge(true);
					//				count++;
					//			}
					//		}
					//	}
					//}
			}
		}
	}
	//delete conv; //Remove excess memory used by convoluter
	//cout << "FOUND " << count << " ridgecells." << endl;
	/*cout << " RESULTS OF CONVOLUTION " << "\n" << endl;
	cout << "--Dx--" << endl;
	cout << outDx << "\n" << endl;

	cout << "--Dy--" << endl;
	cout << outDy << "\n" << endl;

	cout << "--DxDx--" << endl;
	cout << outDxDx << "\n" << endl;

	cout << "--DxDy--" << endl;
	cout << outDxDy << "\n" << endl;

	cout << "--DyDy--" << endl;
	cout << outDyDy << "\n" << endl;*/
	//March through slice to find ridgepoints
	this->grid_points2d = this->m.march2d(this->slice);
	this->grid_lines2d = this->m.constructRidges(this->grid_points2d);

	

}

queue<GridElement<Edge>> RidgeDetector::getRidgePoints2d(){
	return this->grid_points2d;
}

queue<GridElement<Side>> RidgeDetector::getRidgePoints3d(){
	return this->grid_points3d;
}

queue<GridElement<Edge>> RidgeDetector::getRidgeLines2d(){
	return this->grid_lines2d;
}

queue<GridElement<Side>> RidgeDetector::getRidgeLines3d(){
	return this->grid_lines3d;
}

queue<GridElement<Side>> RidgeDetector::getFaceLines(){
	return this->m.temp_lines;
}

RidgeDetector::~RidgeDetector(void)
{
}
